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Abstract 

We calculate quarkonium binding energies using a realistic complex-valued potential for both an 
isotropic and anisotropic quark-gluon plasma. We determine the disassociation temperatures of 
the ground and first excited states considering both the real and imaginary parts of the binding 
energy. We show that the effect of momentum-space anisotropy is smaller on the imaginary part 
of the binding energy than on the real part of the binding energy. In the case that one assumes 
an isotropic plasma, we find disassociation temperatures for the J/^, T and Xb of 1.6 Tc, 2.8 Tc, 
and 1.5 Tc, respectively. We find that a finite oblate momentum-space anisotropy increases the 
disassociation temperature for all states considered and results in a splitting of the p-wave states 
associated with the Xb hi'st excited state of bottomonium. 
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I. INTRODUCTION 



The behavior of nuclear matter at extreme temperatures is now being studied with the 
highest coUision energies ever achieved using the Large Hadron Colhder (LHC) at CERN. 
The ultrarelativistic heavy ion colhsions being studied there will eventually have a center of 
mass energy of 5.5 TeV per nucleon, which is 27.5 times higher than the 200 GeV per nucleon 
energy achieved at the Relativistic Heavy Ion Collider (RHIC) at Brookhaven National 
Laboratory. At RHIC, observations indicated that initial temperatures on the order of twice 
the critical temperature for the quark-gluon plasma phase transition were generated. This 
corresponds to Tq ~ 360 MeV. Assuming that the initial temperature scales with the fourth 
root of the collision energy as predicted by dimensional analysis, one predicts that initial 
temperatures on the order of Tq ~ 4.6 Tc ~ 830 MeV will be generated at the LHC. At such 
high temperatures, one expects to generate a quark-gluon plasma in which the formation 
of quark bound states is suppressed in favor of a state of matter consisting of a deconfined 
plasma of quarks and gluons. 

Suppression of quark bound states follows from the fact that in the quark-gluon plasma 
one expects color charge to be Debye screened ^ . This effect led to early proposals to 
use heavy quarkonium production to measure the temperature of the quark-gluon plasma. 
Heavy quarkonium has received the most theoretical attention since heavy quark states are 
dominated by short rather than long distance physics at low temperatures and can be treated 
using heavy quark effective theory. Based on such effective theories of Quantum Chromo- 
dynamics (QCD) with weak coupling at short distances, non-relativistic quarkonium states 
can be reliably described. Their binding energies are much smaller than the quark mass 
mg ^ Aqcd {Q = c, fe), and their sizes are much larger than l/mq. At zero temperature, 
since the velocity of the quarks in the bound state is small, v <^ c, quarkonium can be 
understood in terms of non-relativistic potential models [3| such as the Cornell potential 
0]. Such potential models can be derived directly from QCD as an effective field theory 
(potential non-relativistic QCD - pNRQCD) by integrating out modes above the scales tuq 
and then mqv, respectively jsf. 

As mentioned above, at high temperature the deconfined phase of QCD exhibits screen- 
ing of static color-electric fields. It is expected that this screening leads to the dissocia- 
tion of quarkonium states, which can serve as a signal for the formation of a deconfined 
quark-gluon plasma in heavy ion collisions [6|. Inspired by the success at zero temperature, 
potential model descriptions have also been applied to understand quarkonium properties 
at finite temperature. The pioneering paper of Matsui and Satz j6| was followed by the 
work of Karsch, Mehr, and Satz 0], which presented the first quantitative calculation of 
quarkonium properties at high temperature. In recent work, more involved calculations of 
quarkonium spectral functions and meson current correlators obtained from potential mod- 
els have been performed js-isll. The results have been compared to first-principle QCD 
calculations performed numerically on lattices 16N24| which rely on the maximum entropy 
method 



A summary and review of the current understanding of potential models is presented 

in 



15[, and different aspects of quarkonium in collider experiments can be found in [28|, 



29| . In recent years, the imaginary part of the potential due to Landau damping has been 



calculated |30H32|. Also, the derivation of potential models from QCD via effective field 



theory methods has been extended to finite temperature (33|. All of the aforementioned 



calculations, however, were performed with the assumption of an isotropic thermal medium. 
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In the last few years there has been an interest in the effect of plasma momentum-space 



anisotropics on quarkonium binding energies for both ground and excited states |34|-|39| . The 



interest stems from the fact that at early times a viscous quark-gluon plasma can have large 
momentum-space anisotropics 40 -4^]. Depending on the magnitude of the shear viscosity, 
these momentum-space anisotropics can persist for a long time (~ 1 - 10 fm/c). The first 
paper to consider the quarkonium potential in a momentum-space anisotropic plasma 34 
considered only the real part of the potential; however, two recent works have extended the 



calculation to include the imaginary part of the potential |36l . l37l |. In this paper we use 



the imaginary part of the potential derived in [37[ and include it in a phenomenological 
model of the heavy quarkonium potential. We then numerically solve the three-dimensional 
Schrodinger equation to find the real and imaginary parts of the binding energies and full 
quantum wavefunctions of the charmonium and bottomonium ground states, as well as the 
first excited state of bottomonium. We present data as a function of the temperature and 
are able to identify the full effect of the isotropic and anisotropic potentials on these states. 

We compare our results with a recent analytic estimate of the imaginary part of the 
binding energy by Dumitru [49[. We show that, for an isotropic plasma, the imaginary 
part of the binding energy is approximately linear in the temperature for temperatures 



near the phase transition in agreement with Ref. j49|. However, in the case of the J/if)^ 
we find a significantly smaller slope for the imaginary part of the binding energy as a 
function of temperature than predicted by Ref. 4^. The discrepancy most likely arises 



from the fact that Ref. [49[ assumed Coulombic wavefunctions. The potential used here 
includes modifications at both intermediate and long ranges, which causes the numerical 
wavefunctions to not be well approximated by Coulombic wavefunctions. In addition, our 
wavefunctions are complex with the imaginary part growing in magnitude as the temperature 
is increased. This effect was ignored by the assumption of Coulombic wavefunctions in 
Ref. [49|. We find that when the states are small and dominated by the screened Coulomb 
potential, the imaginary part of the binding energy increases approximately linearly with 
the temperature; however, as the size of the bound state increases, the scale set by the string 
tension dominates and the imaginary part of the binding energy increases more slowly with 
increasing temperature. 

The structure of this paper is as follows: In Sec. [ITl we review the potential introduced 



in Ref. [35| and extend it to include the imaginary part of the potential derived in Ref. [37 
In Sec. we review the numerical method that we use to solve the three-dimensional 
Schrodinger equation. In Sec. llVt we present our numerical results for the real and imaginary 
parts of the binding energies of the charmonium and bottomonium ground states and first 
excited state of bottomonium. In Sec. |Vl we state our conclusions and give an outlook for 
future work. Finally, in an appendix we present numerical benchmarks and tests of the code 
used here in order to demonstrate its convergence and applicability to the problem at hand. 



II. SETUP AND MODEL POTENTIAL 



In this section we specify the potential we use in this work. We consider the general 
case of a quark-gluon plasma which is anisotropic in momentum space. In the limit that 
the plasma is assumed to be isotropic, the real part of the potential used here reduces to 
the model originally introduced by Karsch, Mehr, and Satz (KMS) 0] with an additional 
entropy contribution [3^ and the imaginary part reduces to the result originally obtained 
by Laine et al [l^l- To begin the discussion we first introduce our ansatz for the one-particle 
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distribution function subject to a momentum-space anisotropy. 



A. The anisotropic plasma 

The phase-space distribution of gluons in the local rest frame is assumed to be given by 



the following ansatz [3J, |50|-|53 



/(x, p) = /i3o [Vp' + e(p ■ n) VpLd) , (1) 

where Phard is a scale which specifies the typical momentum of the particles in the plasma and 
can be identified with the temperature in the limit that ^ = 0. Thus, /(x, p) is obtained from 
an isotropic distribution /iso(|p|) by removing particles with a large momentum component 
along n, the direction of anisotropy. In this paper, we will restrict our consideration to a 
plasma that is close to equilibrium. This is motivated by the fact that in a heavy-ion collision, 
quarkonium states are expected to form when the temperature has dropped to (1-2) Tc. At 
such temperatures the plasma may have partly equilibrated/isotropized. Additionally, this 
means that we can assume that the function /iso(|p|) is a thermal distribution function. 
The parameter determines the degree of anisotropy, 

where = p ■ n and p^ = p — n(p ■ n) denote the particle momentum along and perpen- 
dicular to the direction n of anisotropy, respectively. If ^ is small, then it is also related to the 
shear viscosity of the plasma. For example, for one-dimensional boost-invariant expansion 



governed by Navier-Stokes evolution j45|, |47|, |48|, |5J] one finds 



i-p, (3) 

ITS 

where T is the temperature, r is the proper time (and 1/r is the Hubble expansion rate), and 
?7/s is the ratio of shear viscosity to entropy density. In an expanding system, non- vanishing 
viscosity (finite momentum relaxation rate) implies an anisotropy of the particle momenta 
which increases with the expansion rate 1/r. For rj/s ~ 0.1 - 0.2 and tT ~ 1 - 3 one finds 
that ^ — 1. In general, one can relate ^ to the longitudinal and transverse pressures in the 
plasma and it is possible to derive dynamical differential equations which govern its time 



evolution similar to viscous hydrodynamics |47|, |48 

We point out that in this paper we restrict ourselves to solving the time-independent 
Schrodinger equation, i.e. we assume that the plasma is at a constant hard momentum scale 
Phard and anisotropy C,. This approximation is useful if the time scale associated with the 
bound state, ~ l/l-Ebmdl, is short compared to the time scales over which Phard and ^ vary. 
Indeed, for sufficiently large quark mass mg this condition should be satisfied. 

B. The model potential 

Lacking knowledge of the exact heavy-quark potential at finite temperature, different 
phenomenological potentials and lattice-QCD based potentials have been used to study 
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quarkonium binding energies in the quark-gluon plasma. To start, we decompose the po- 
tential into real and imaginary parts, V = Vn + iVi. The model for the real part of the 
potential we use was obtained in Ref. j35|. The analytic calculation of the imaginary part 



was performed in Refs. [30|, l37|, |55|. The real part is given by 



Va(r) = -- (1 +/ir)exp(-/ir) + — [1 



where 



exp(— /ir)] —ar exp(— /ir) 



3 + cos 29 



0.8 a 

2 ' 



(4) 



(5) 



m£) 16 

with rriD = (1.4)^ ■ Nc{l + Nf/6) Anaspl^^^/S being the isotropic leading-order Debye mass 



adjusted by a factor of (1.4)^ to take into account higher-order corrections [56|. The coupling 
a folds in a factor of Cp = {N^ — l)/{2Nc), i.e. a = Cpcts, where as = g1/{AT{) is the 
canonically defined strong coupling constant. We have taken A'^c = 3 and assumed Nf = 2 
which is appropriate for the temperature range considered herein. The first term in (j4]) is 
a screened Coulomb potential with an entropy addition. The second and third terms are 
a screened linear potential associated with confinement in the low temperature limit. The 
last term in is a relativistic correction which is critical for obtaining accurate binding 
energies in the low temperature limit. For the string tension, we fix a = 0.223 GeV and for 
the strong coupling constant we fix a = 0.385.^ 
The imaginary part is given by ^ 



Vi{r) 



-aT 



where r = m£,r and 



dz- 



z 



Mr,0) = I dz 

'0 





oo 



Z^ + 1)2 

Z 



1 - 



simzr] 



1-5 

[z' + l)''\ 2 



zr 



sin^^^^ , n 



dz- 







1-3 



z r 



cosH 



l-3cos^e)G{r,z) 



zr 



[^2 + 1)3 

with 6 being the angle from the beam direction and 

rz cos(fz) — sm(fz) 
G{r,z) = — — 



'l-3cos^9)G{f,z) 



(6) 

(7) 
(8) 

(9) 



(10) 



The short range part of V r is based on a leading order hard-loop perturbation theory 
calculation presented in Ref. [3J|. Vi is also obtained from a leading order perturbative 
calculation [37|. Being a leading order calculation one may wonder about higher order 



Since as runs logarithmically and therefore has small variation in the temperature ranges shown, we will 
ignore the running of the coupling here. Incorporating this effect would be straightforward, however, a 
model of the behavior of at large scales would be required in order fit zero temperature properties of 
the states considered here. 
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corrections. One expects that the leading order calculation pQCD would receive large cor- 
rections at low temperatures (T < 10 Tc) since the running coupling becomes large {qs > 1). 
For the coupling used above as = 0.29 one finds Qs = y/A¥a^ =1.9. This means that the nor- 
mal scale hierarchy, QgT < T, implicit in the hard-loop resummation becomes inverted.^ We 
therefore need to supplement the leading order pQCD calculation with a non-perturbative 
contribution. For the real part we do this by including a long-range screened linear contri- 
bution that is modified to include an entropy contribution 3J]. In the isotropic limit the 
resulting form of the real part potential is in good agreement with lattice data for the heavy 
quark potential [56[. For the imaginary part we currently do not have non-perturbative 
input from lattice calculations with which to constrain the long range part; however, we 
note that calculations of the real and imaginary parts of the potential using the AdS / CFT 
correspondence to calculate the corresponding potential in large t' Hooft coupling limit of 
A/" = 4 Supersymmetric Yang- Mills yield similar results to those obtained using perturbative 
6l|. For more information about the relevant scales and limitations of the current 



QCD [381, 

approach we refer the reader to Sec. Ill of Ref. 34 . 

Regarding the length scales which are relevant, we note that the short range part of the 
potential is appropriate for describing wavefunctions which have l/(r) < 0{mD) while the 
long range part is relevant if l/(r) > O^mo)- Using the form of the real potential listed 
above, one finds that the distance scale at which medium effects become large is roughly 
given by r > rmed ~ Tc/ (2T) fm corresponding to rmed ~ 0.25 fm at 2 Tc jl^l • Numerically, 
the isotropic Debye mass used herein is m^) ~ 3phard, corresponding to rriD ~ 1.2 GeV at 
Phard = 2Tc. As showu in Ref. jl^l Fig. 4, using the real part of the potential listed above, the 
RMS radius of the J/\l/ state is approximately 0.8 fm at 2 Tc corresponding to 1/ (r) ~ 250 
MeV, which makes the screening of the long range part of the potential crucially important 



for fixing the binding energy in this case. For the case of the T one sees also from Ref. [34 
Fig. 4 that the RMS radius of the T is approximately 0.25 fm corresponding to 1/ (r) ~ 800 
MeV. We note importantly that for the T, due to its relatively small size, the bulk of the 
medium effect comes from the temperature dependence of lim^-i-oo V = Voo (see Fig. 3 of 
Ref. 3J]). In closing, one finds that for both the J/\l/ and T that correct modeling of both 
the short and long range parts of the potential are critical for obtaining the temperature 
dependence of these states. As mentioned above, here we extend the results in [3j| to include 
the imaginary part of the potential. We note that one finds that RMS radii of the states are 
only weakly affected by inclusion of the imaginary part of the potential, allowing us to use 
the estimates above as a rough guide for understanding the relevant scales. 



C. Analytic estimate in isotropic case 



In a recent paper [49 



Dumitru made an estimate of the effect of the imaginary part of 
the potential on the imaginary part of the binding energy of a quarkonium state. For this 
estimate Dumitru assumed a Coulomb wavefunction for the quarkonium state and computed 
the expectation value of the imaginary part of the potential exactly in the case of an isotropic 



^ We note that for temperatures T > 2Tc NNLO perturbative calculatfons of QCD thermodynamics based 
on hard-thermal-loop resummation of QCD agree quite well with available lattiee data even though is 



large [57h60| 
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plasma. The result obtained was 



a ttIq (1 — K^)-^ a niQ 

When plotted in the temperature range between Tc and 3Tc the result above is approximately 



linear for both the J/i/j and T [49|]. For charmonium with mq = 1.3 GeV and using the 
values given for a and m^j in the previous subsection, we obtain a slope consistent with 
r oc (0.08 GeV) T/T^ at T = 0.3 GeV. Similarly, for bottomonium with mq = 4.7 GeV we 
obtain a slope consistent with F oc (0.05 GeV) T /Tc. We note these here for later comparison 
with numerical results presented in the results section. 



III. NUMERICAL METHOD 

To determine the wavefunctions of bound quarkonium states, we solve the Schrodinger 
equation 

i^'0„(x) = i?„0„(x) , 

H = + r(x)+mi+m2, (12) 

Irrifi 

on a three-dimensional lattice in coordinate space with the potential given hj V = Vr^ + iVi 
where the real and imaginary parts are specified in Eqs. (jl]) and ([6]), respectively. Here, 
mi and m2 are the masses of the two heavy quarks and m/j is the reduced mass, mji = 
17111712/ {mi +777.2). The index v on the eigenfunctions, 0„, and energies, E^, represents a 
list of all relevant quantum numbers, such as 77, /, and 777 for a radial Coloumb potential. 
Due to the anisotropic screening scale, the wavefunctions are no longer radially symmetric 
if ^ 7^ 0. Since we consider only small anisotropies we nevertheless label the states as IS* 
(ground state) and IP (first excited state), respectively. 

To find solutions to Eq. f[T^ . we use the finite difference time domain method (FDTD) [62|, 



63|. In this method we start with the time-dependent Schrodinger equation 



i^7/.(x,t) = ^7/;(x,t), (13) 
which can be solved by expanding in terms of the eigenfunctions, 0„: 

^{^,t) = J2cvM^)e-'''-' . (14) 

V 

If one is only interested in the lowest energy states (ground state and first few excited states) 
an efficient way to proceed is to transform (fT3l) and (fT^ to Euclidean time using a Wick 
rotation, r = it: 

^7A(x,r) = -if7A(x,r), (15) 

and 

V.(x,r) = 5^c„</.„(x)e-^^^. (16) 



For details of the discretizations used etc. we refer the reader to Refs. 62, 63 
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A. Finding the ground state 



By definition, the ground state is the state with the lowest energy eigenvalue, Eq. There- 
fore, at late imaginary time the sum over eigenf unctions f|T6|) is dominated by the ground 
state eigenfunction 



Due to this, one can obtain the ground state wavefunction, 0O; and energy, Eq, by solving 
Eq. (fT5|) starting from a random three-dimensional wavefunction, -i/^initiaKx, 0), and evolving 
forward in imaginary time. This initial wavefunction should have a nonzero overlap with 
all eigenfunctions of the Hamiltonian; however, due to the damping of higher-energy eigen- 
functions at sufficiently late imaginary times we are left with only the ground state, 0o(x)- 
Once the ground state wavefunction (or any other wavefunction) is found, we can compute 
its energy eigenvalue via 



To obtain the binding energy of a state, -Ew^bind; we subtract the quark masses and the 
real part of the potential at infinity 



For the isotropic KMS potential the last term is independent of the quantum numbers v and 
equal to 2(j/m£). In the anisotropic case, however, this is no longer true since the operator 
Voq{9) carries angular dependence, as discussed above. Its expectation value is, of course, 
independent of 9 but does depend on the anisotropy parameter ^. 

B. Finding the excited states 

The basic method for finding excited states is to first evolve the initially random wave- 
function to large imaginary times, find the ground state wavefunction, 0o, and then project 
this state out from the initial wavefunction and re-evolve the partial-differential equation in 
imaginary time. However, there are (at least) two more efficient ways to accomplish this. 
The first is to record snapshots of the 3d wavefunction at a specified interval Tsnapshot during 
a single evolution in r. After having obtained the ground state wavefunction, one can go 
back and extract the excited states by projecting out the ground state wavefunction from 
the recorded snapshots of r). 

An alternative way to select different excited states is to impose a symmetry condition 
on the initially random wavefunction which cannot be broken by the Hamiltonian evolution. 
For example, one can select the first excited state of the (anisotropic) potential by anti- 
symmetrizing the initial wavefunction around either the x, y, or z axes. In the anisotropic 
case this trick can be used to separate the different polarizations of the first excited state of 
the quarkonium system and to determine their energy eigenvalues with high precision. This 
high precision allows one to more accurately determine the splitting between polarization 
states which are otherwise degenerate in the isotropic Debye-Coulomb potential. 




(17) 




(18) 




(19) 
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FIG. 1: Real and imaginary parts of the charmonium ground state { J/ip) binding energy as a 
function of phard- Both isotropic ^ = and anisotropic ^ = 1 cases are shown. The left panel shows 
full temperature range and the right panel focuses on the region where the real and imaginary 
parts become comparable. See text for parameters such as lattice size, lattice spacing, etc. 



Whichever method is used, once the wavefunction of an excited state has been determined 
one can again use the general formulas (|T8l) and ( |T9|) to determine the excited state binding 
energy. For code benchmarks and tests see App. \M 



IV. RESULTS AND DISCUSSION 

In this section we present results for isotropic = 0) and anisotropic = 1) binding en- 
ergies for charmonium (J/ip), bottomonium (T), and the first excited state of bottomonium 
(Xb) as a function of the hard momentum scale phard- We will first assume that phard is held 
constant and vary the anisotropy parameter. Note increasing ^ results in a decrease in the 
density since n oc Phard/ Vl + ^ ^M- This reduced density results in less Debye screening and 
thus a more strongly bound state. We therefore expect that states with large anisotropy will 
have increased binding energies compared to the isotropic states. One could imagine holding 
another thermodynamic property such as the number density or energy density constant as 
one changes the anisotropy parameter. We will return to this issue at the end of this section 
and show that the results in these cases can be obtained from a simple rescaling of the results 
presented below. In all plots shown, we assume = 192 MeV and fix the imaginary-time 
step in the numerical algorithm to be At = a^/8 where a is the spatial lattice spacing. 



A. Results as a function of the hard momentum scale 

In Fig. [1] we plot the binding energy of the charmonium ground state {J/ip) as a function 
of Phard- ^OT this figurc, wc uscd a lattice size of 256'^ with lattice dimension of L = 25.6 
GeV~^ and a lattice spacing of a = 0.1 GeV~^. For the charmonium mass, we used rUc = 1.3 
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FIG. 2: Real and imaginary parts of the bottomonium ground state (T) binding energy as a 
function of Phard- Both isotropic ^ = and anisotropic ^ = 1 cases are shown. See text for 
parameters such as lattice size, lattice spacing, etc. 



GeV. In Fig. [T] we show both the real part (black line with filled circles) and the imaginary 
part (red line with filled squares) of the isotropic ground state binding energy. Comparing 
these two curves, we see that the imaginary part of the binding energy becomes comparable 
to the real part at phard ~ 1.63 Tc. In contrast, in the anisotropic case (^ = 1) we find that 
the intersection between the imaginary (blue line with open triangles) and real parts (green 
line with open diamonds) occurs at Phard ~ 1.88 T^. In the range between 1 and 3 Tc we 
obtain a slope of 4.9 x 10~^ GeV for the imaginary part of the binding energy when ^ = 
and 6.4 x 10~^ GeV when ^ = 1. In the isotropic case Dumitru's perturbative calculation 
(49I gives a slope of 8 x 10"^ GeV. Our method is non-perturbative since we don't assume 
perturbations around Coulomb wave functions, so one should not be surprised to see some 
important differences. 

In Fig.[2]we plot the binding energy of the bottomonium ground state (T) as a function of 
Phard- For this figure, we used a lattice size of 256^ with lattice dimension of L = 25.6 GeV~^ 
and a lattice spacing of a = 0.1 GeV~^. For the bottomonium mass, we used rrib = 4.7 GeV. 
In Fig. [2] we show both the real part (black line with filled circles) and the imaginary part 
(red line with filled squares) of the isotropic ground state binding energy. When ^ = 0, we 
see that the imaginary part of the binding energy becomes comparable to the real part at 
Phard ~ 2.8 Tg. In the anisotropic case = 1) we find that the intersection between the 
imaginary (blue line with open triangles) and real parts (green line with open diamonds) 
occurs at approximately 3.5 Tc. For ^ = 0, in the range between 1 and 4 Tc we obtain a 
slope of 2.8 X 10"^ GeV for the imaginary part of the binding energy. In the anisotropic 
case (^ = 1) we find a slope of 4.2 x 10~^ GeV. We can once again compare to the analytic 
result of Dumitru j49| which gives an isotropic slope of 5 x 10"^ for the T. Once again, the 
numbers are roughly in agreement. 

In Fig. [3] we plot the binding energy of the first p-wave excited state of bottomonium (xb) 
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FIG. 3: Real and imaginary parts of Xb binding energy as a function of J'hard- Both isotropic ^ = 
and anisotropic ^ = 1 cases are shown. See text for parameters such as lattice size, lattice spacing, 
etc. 



as a function of Phard- For this figure we used a lattice size of 256'^ with lattice dimension of 
L = 38.4 GeV~^ and a lattice spacing of a = 0.15 GeV~^. For the bottomonium mass, we 
used rrib = 4.7 GeV. As was the case with the bottomonium ground state we see an increase 
in the real part of the binding energy with increasing anisotropy. Most importantly, we 
find that there is an approximately 60 MeV splitting between the = and Lz = ±1 
states with the states with = ±1 having the lower binding energy. We would therefore 
expect fewer = ±1 states of the xt to be produced in an anisotropic plasma. Determining 
precisely how many fewer would be produced requires knowledge of the time evolution of 
the momentum scale Phard and anisotropy ^. 



B. Fixing number density or energy density 

As mentioned in the beginning of this section, when one is working in a non-equilibrium 
setting it is necessary to specify which quantities are held fixed. In equilibrium, it is sufficient 
to specify the temperature. The temperature then uniquely determines the number density, 
energy density, etc. In the previous subsection we presented results obtained when one holds 
the hard momentum scale Phard fixed while varying the anisotropy parameter ^. Doing so, 
however, results in different number densities and energy densities for different anisotropies 
(^). Here we discuss how to fix either the number or energy density by adjusting Phard 
appropriately. We first demonstrate this in the case of the number density and show that 
for small anisotropy the scalings required to fix the number density or energy density are 
practically identical. We then present results for the binding energies of the states we are 
interested in for the case of fixed number density, since in this paper we concentrate on 
anisotropies which are small enough that the difference between the cases of fixed number 
density and fixed energy density is numerically very small. 
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FIG. 4: Real and imaginary parts of the charmonium ground state (J/^) binding energy as a 
function of temperature assuming fixed number density. Both isotropic ^ = and anisotropic 
^ = 1 cases are shown. See text for parameters such as lattice size, lattice spacing, etc. 



The number density as a function of ^ and Phard can be calculated for an arbitrary isotropic 
distribution /iso j52| 

I i- \ ''^isolPhardJ /on\ 

n(e,Phard) = ^JYTl ' ^ ^ 

where n\^o is the number density associated with the isotropic distribution function /iso via 

^^^/iso(|p|,Phard) • (21) 

Since riiso only contains one dimensionful scale, by dimensional analysis we have riiso oc p^ard- 
In order to keep the number density (l20l) fixed as one changes ^, one can adjust Phard by 
requiring 

Phard = (1 + 0^^^ ^ [fixed number density] , (22) 

where T is the corresponding isotropic scale (temperature) which gives the target number 
density when = 0, i.e. niso(7')- 

Similarly, the energy density as a function of ^ and Phard can be calculated for an arbitrary 
isotropic distribution /iso js^] 



^ (^, Phard) = ^(O^iso(Phard) , (23) 

where £^iso is the energy density associated with the isotropic distribution function /iso and 

7^(0 = -fA: + ^^^^V (24) 



2Vl + e 

Since £^iso only contains one dimensionful scale, by dimensional analysis we have £^iso cx: p^ard 
and we can fix the energy density to the corresponding isotropic energy density with scale 



12 




FIG. 5: Real and imaginary parts of bottomonium ground state (T) binding energy as a function 
of temperature assuming fixed number density. Both isotropic ^ = and anisotropic ^ = 1 cases 
are shown. See text for parameters such as lattice size, lattice spacing, etc. 



T by requiring 

Phard = T/[7^(0]'/' [fixed energy density] . (25) 



The scalings for fixed number density (12211 and fixed energy density ( 1251) are different; 
however, in the limit of small anisotropies the scalings are very close. Expanding to quadratic 
order, one finds 

^ = 1 + 1^- + Oi^) [fixed number density] , (26a) 

^ = l + l^-^e + Oie) [fixed energy density] , (26b) 

which agree at linear order and differ by 7.4% in the quadratic coefficient. One finds that, 
when including all orders in the expansion, the right hand sides of ( 12 6 a!) and ( I26bl) differ 



by only 0.25% at ^ = 1. Therefore, for the range of anisotropies considered here, the two 
scalings are functionally equivalent. We will therefore only present results for fixed number 
density with the understanding that the fixed energy density results are indistinguishable 
by the human eye. 

In Figs, m [5l and O we show the binding energies which result from the fixed number 
density rescaling of the horizontal axes of Figs. (H [2l andO As can be seen from Figs. IH |5l and 
El requiring fixed number/energy density weakens the effect of anisotropies on the ground 
state binding energies. In the case of the ground states of charmonium and bottomonium 
shown in Figs. H] and E] we find that the splitting between the ^ = and ^ = 1 cases at the 
critical temperature is approximately 50 MeV in both cases. 

Finally, we emphasize that in the case of the first excited states of bottomonium shown 
in Fig. [6] the splitting between the = and Lz = ±1 states is unaffected by the rescaling 
since we have ^ = 1 for both states. Therefore, one has a relatively clean observable that 
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FIG. 6: Real and imaginary parts of the Xb binding energy as a function of temperature assuming 
fixed number density. Both isotropic = and anisotropic ^ = 1 cases are shown. See text for 
parameters such as lattice size, lattice spacing, etc. 



is sensitive to plasma anisotropies regardless of the quantity which is assumed to be held 
fixed. 



V. CONCLUSIONS 

In this paper we have presented first results on the effect of including both the real 
and imaginary parts of the heavy quarkonium potential on the binding energies of the 
charmonium ground state {J/ip), the bottomonium ground state (T), and the first p-wave 
excited state of bottomonium (xb)- We did this by numerically solving the three-dimensional 
Schrodinger equation for the complex potential given by Eqs. (jlj) and ([6]). This enabled us 
to extract both the real and imaginary parts of the binding energies for the states. Using 
our model potential, we investigated both isotropic and weakly anisotropic plasmas. We 
found that, there can be a sizable effect of momentum-space anisotropy on both the real 
and imaginary parts of the quarkonium binding energy. One can estimate the disassociation 
temperature of the states by determining the temperature at which the real and imaginary 
parts of the binding energy become the same. Using this criteria, in the isotropic case we 
estimate the J/ip, T and Xb to have disassociation temperatures of 1.6 T^, 2.8 T^, and 1.5 
Tc, respectively. We note, however, that even prior to these disassociation temperatures the 
states will be suppressed due to the exponential decay of the states with a rate related to the 
imaginary part of the binding energy. We plan to investigate the phenomenological impact 
of our results on the time evolution of quarkonium decay in a future publication. 

In the case of a plasma with a finite momentum-space anisotropy, we presented results for 
both fixed hard momentum scale and fixed number density. Our results demonstrate that the 
corresponding anisotropic states have a higher binding energy in accordance with pre vious 



results that employed only the real part of the quarkonium potential used herein [35|. We 
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showed that, for small anisotropy, fixing the number density and fixing the energy density 
gives results which are the same to within less than a fraction of a percent. We demonstrated 
that fixing the number density reduces the effect of anisotropy compared to the case of fixing 
the hard momentum scale, but does not completely remove the effect of momentum-space 
anisotropy on the binding energies. Finally, we emphasized the importance of the finite- 
anisotropy splitting between the Xb states with = and = ±1. This splitting is 
independent of whether one fixes the hard momentum scale, number density, or energy 
density. Therefore, this splitting represents a possible observable which could be used to 
determine the time-averaged plasma anisotropy parameter. 

Looking forward, to fully assess the phenomenological impact of plasma momentum-space 
anisotropies on quarkonium states requires the convolution of the results presented here with 
the space-time evolution of the hard momentum scale and anisotropy parameter. A method 
for determining the dynamical evolution of these parameters has recently been determined 



47|, |48|. In addition, since these works show that ^ can become large, it will be necessary to 
investigate the effect of large anisotropies on quarkonium binding energies. The calculations 
necessary to address these questions are currently underway. 
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Note added 

In arXiv versions 1-3 there was a mistake in the final results for the imaginary part of 
the binding energies. The mistake stems from the fact that we had subtracted the full 
complex- valued potential at infinity Voo] however, formally only the real part of Voo should 
be subtracted since the imaginary part of Voo is related to heavy quark damping in the 
plasma which is physically relevant. As a consequence all imaginary parts of the binding 
energies are changed and we have updated all figures. While the results are qualitatively 
similar, the key change is that the imaginary part of the binding energy now has a stronger 
dependence on the anisotropy parameter, ^, in most cases. 



Appendix A: Numerical Tests 
1. Convergence Test 

In this appendix we present some convergence data for a particular state in order to 
demonstrate the approach to the continuum limit. In Fig. [7] we show both the real and 
imaginary parts of the bottomonium ground state binding energy for three different lattice 
spaces of 128"^, 256^, and 512"^. For each of the runs, the lattice size was fixed to L = 25.6 
GeV~^ with the lattice spacing in each case given by a = 0.2, a = 0.1, and a = 0.05 GeV~^, 
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Imaginary Part 




FIG. 7: Real and imaginary parts of the bottomonium binding energy for three different lattice 
sizes of 128^2563, and 512^. 



respectively. We chose ^ = 1, rrib = 4.7 GeV and used an imaginary-time step given by 
At = in each case. 

As can be seen from this figure, there is a larger effect due to reducing the lattice space 
on the real part than the imaginary part. This is to be expected since the real part contains 
a divergence at the origin, whereas the imaginary part is regular. In the case of the real 
part at 3 Tc, we see an approximately 8.0% correction when going from 128'^ to 256^ and a 
2.6% correction when going from 256^ to 512^. The corrections in the case of the imaginary 
part are 2.4% and 0.16%, respectively. Therefore we see that the 256^ runs presented in the 
body of the text are reliable up to corrections on the order of 3%. 



2. Harmonic Oscillator with complex spring constant 

In this part of the appendix we explore the ability of the FDTD algorithm to handle 
the case of complex potentials. We investigate a simple one dimensional test case which 
consists of solving the Schrodinger equation for a particle of mass m in a harmonic oscillator 
potential which has a complex spring constant k. We first derive the analytic solution and 
then compare the output of the code and the analytic solution. 

Our goal is to solve the time-independent Schrodinger equation for a particle of mass m 
which is bounded by the quadratic potential V{x) = kx^/2, in the case that k is complex. 
Before proceeding, we review the solution in the case that k is real, writing it first in 
terms of Parabolic Cylinder Functions and then showing how these reduce to the Hermite 
Polynomials. 
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a. Review of the case of a real spring constant 

We are interested in solving the time independent ID Schrodinger equation for the har- 
monic oscillator potential 

To begin with, we introduce the variables uj = -y/^, a = (^)^''^, = ^J^d u = ax 
which allow us to write (lAlll compactly as 

"^'"^ {u^-X)i> = 0. (A2) 



b. Asymptotic Behavior 



We now want to find solutions for our wavefunction at |n| — )■ oo. In the limit |n| >> 1 
Eq. flA2p becomes 

hm P^=u'^P, (A3) 
which has an asymptotic solutions of the form 

lim 4j{u) = yl^^Pe-"'/^ (A4) 

or 

lim ip{u) = 5u'?e"'/^ (A5) 

\u\—^oo 

with A and B being arbitrary constants. Since ip must remain finite as \u\ oo, this requires 
that we discard the second solution. Requiring that the wavefunction be single valued for 
negative u implies that p must be an integer.^ This yields an asymptotic solution of the 
form 

lim ij{u) = AuPe-'''/^ (A6) 

\u\—^oo 

where p is integer-valued. 



c. Solution in terms of Parabolic Cylinder Functions 

We now define new variables a = —A/2 and b = ^/2u such that u'^ = 6^/2 and 2cP/dlP' = 
d'^/du'^. Substituting these into Eq. (IA2|) gives 

^y + aV = 0. (A7) 



db^ \A 

The solution to this differential equation is given by the Parabolic Cylinder Function D from 



Chapter 19 of Abramowitz and Stegun [65|] ilj{a,b) = U{a,b) = -D„a„if,(&). If U{a,b) is a 



^ Note that one could move the cut along the negative u axis into the complex plane, so the more properly- 
stated requirement is that the wavefunction be single- valued everywhere in the complex plane. 
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solution, then so are U{a, —b),U{—a,ib), and U{—a, —ib). This leaves us with the general 
solution 



iP{u) = AU (^~,V2u^ +BU (^^,iV2uy (AJ 



d. Matching Asymptotic Behavior 



Now we analyze the asymptotic behavior of U when b ^ a. For 6 ^ a, we have (65 

f/(a, b) ~ e-3^'r"-^ (1 - C(r2)). (A9) 

We see that in the limit 6 — )■ oo that U (a, ife) approaches +oo. Since the wavefunction should 
be finite at — >■ oo this requires 5 = in Eq. (lASp . We can solve for A by matching to the 
asymptotic form given in Eq. ( ]A6P 



A 1 

-2 + 2 



-p. (AlO) 



where p must be an integer as discussed previously. Solving for A gives us A = 2p + 1, 
which tells us that A is discrete. This leaves us with the solution to the quantum harmonic 
oscillator differential equation in terms of Parabolic Cylinder D functions 

^{u) = Au(-'^^^,V2u] . (All) 



e. Connection to Hermite Polynomials 



When n is a non-negative integer, U{—n — 1/2, b) is expressible in terms of Hermite 
Polynomials j65| 

U{-n - 1/2, b) = 2~-^^e--^''H^ {^-^ , (A12) 

where Hn{x) is a Hermite polynomial. Using this we find that Eq. (lAlip can be expressed 
as 

iPp{u) = Apc"^'^' Hp{u). (A13) 

This is the standard textbook form of the quantum harmonic oscillator wavefunctions. We 
now extend the general solution flASp to the case where k can be complex. 



3. Extending the solution to the case of a complex spring constant 

We now examine the behavior of our solution when the spring constant k is complex. To 
begin, we note a symmetry of Eq. f lA2p : namely, if we find that a solution for a given k is 
given by then the solution for the complex conjugate spring constant k* will be given by 
if)*. Therefore, it suffices to solve the equation in only half of the complex k plane. In order 
to find the specific solution, we take the general solution in Eq. ( lASP and re-examine its 
asymptotic behavior. Depending on the angle of k in the complex plane, we find that one 
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needs to set either A or 5 in Eq. (lASP to zero. For large values of m >> \a\ we have from 

Ref. lei 



1„2 



[/(a,u) ~ e~3« u~''"2, (A14) 

where we have dropped a factor of a/2 for simplicity. For complex-valued m = /3 + ^7 with 
real-valued /3 and 7, we have three cases that determine the asymptotic convergence of 
U{a,u): /3 > 'J, f3 < 'J, and /5 = 7. 

a. Case A: Real part greater than imaginary part (13 > 

In this case, we must require B = in Eq. ( ]A8I) and we have 

^(m) = [/(a,/3 + i7) ~ e-^('^+*^)'M-'^-5 (A15) 

+2^/37-72 -a- i 



-1/32 -ijfl^ 1^2 



which converges since /3 > 7. 

6. Case B: Real part less than imaginary part (f3 < ^) 

In this case, we require A = in Eq. (1A8I1 and we have 

U{-a,iu) = U{-a,i(3--f) (A16) 

= g-i(-/32-2i/37+7')^a-i 



which converges since < 7. 

c. Case C.' Real part equal to imaginary part (f3 = ^) 

In this case, the solution diverges; however, as we will show below, /3 = 7 corresponds to 
a purely repulsive potential so one could expect such singular behavior. Since here /3 = 7, 
we can rewrite 

u = p + if3= |n|e*"/^ (A17) 
Also, since u = [mk/h^Y^^x and k = |/c|e*^, we have 

u = k^'\m/h^Y'^x (A18) 
= e''^/\m\k\/h^Y/^x (A19) 

Equating the imaginary parts of Equations (1A17I) and (1A19I) . we have 

gW4 ^ ^^20) 

which means 6 = ir and k = —l\k\. This corresponds to a repulsive spring constant which 
has no bound solutions. 
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FIG. 8: Plot showing real and imaginary parts of the wavefunction for a complex harmonic oscillator 
potential with k = 1 + i. The solid lines are the analytic results and the circles are sampled points 
obtained from the numerical solution using the FDTD method [63 1. 



4. Comparison of Numerical Solution and Analytic Solution 

In this section, we will examine two different complex values of k. Since u = ax where 
a = {mk/h?y^^, the solution depends on the real and imaginary parts of the spring constant 
k. This determines whether A = or 5 = in Eq. (lASp as discussed above. Below, we use 
natural units with h = c = 1 and take m = 1 GeV. 



a. Case k = 1 + i 

In this case, we have k = e^'^^^\k\ and u = e™^^^^Ur where Ur = (mlkl/h^Y^^x. Therefore, 
/3 > 7 and we must use case A. The solution becomes 

i:{u) = Au(^-^,V2uy (A21) 

In Fig. |8]we plot the result for the ground state which corresponds to A = 1 with the con- 
stant A fixed to require a normalized wavefunction. The solid lines are the analytic results 
and the circles are sampled points obtained from numerical solution using the FDTD method 



63|.^ As we can see from this figure, the FDTD algorithm is able to obtain very good agree- 
ment in both the real and imaginary parts of the wavefunction. We can also compare the 
ground state energy predicted analytically by Eq. flMT]) which is Eq = 0.549342 + 0.227545i 



with the FDTD algorithm's result. The FDTD algorithm, using 200 points distributed with 
a step size of 0.05 and a convergence tolerance of 10~^, gives Eq = 0.549251 + 0.227479z, 
which represents an accuracy of approximately 0.01%. 



^ The wavefunction can be rotated by an arbitrary complex phase e** without affecting the probability 
amplitudes. In practice, the code converges to a different random phase angle during each run. In order 
to compare to the analytic results for the real and imaginary parts of the wavefunction we manually rotate 
the numerically determined wavefunctions such that Im[i/'(a; = 0)] = 0. 
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FIG. 9: Plot showing real and imaginary parts of the wavefunction for a complex harmonic oscillator 
potential with k = —i. The solid lines are the analytic results and the circles are sampled points 
obtained from numerical solution using the FDTD method 63|. 



b. Case k = —i 

In this case, we have k = e*^'^/^|/c| and u = e^^'^^^^Ur where Ur = {m\k\/h'^Y^^x. Therefore, 
(3 < 'J and we must use case B. The solution becomes 

ip{u) = 5 [/ Q, iV2v^ . (A22) 

In Fig. [9] we plot the result for the ground state which corresponds to A = 1 with B 
fixed to require a normalized wavefunction. The solid lines are the analytic results and the 



circles are sampled points obtained from numerical solution using the FDTD method [63 
As we can see from this figure, the FDTD algorithm is able to obtain very good agreement 
in both the real and imaginary parts of the wavefunction. We can also compare the ground 



state energy predicted analytically by Eq. f lA22p which is Eq = 0.353553 — 0.353553i with 



the FDTD algorithm's result. The FDTD algorithm, using 200 points distributed with step 
size of 0.05 and a convergence tolerance of 10"^, gives Eq = 0.353522 — 0.353478i, which 
represents an accuracy of approximately 5 x 10"'^ %. 

The results of the two cases presented in this appendix show that the FDTD algorithm is 
able to obtain accurate wavefunctions and eigenvalues even in the case that the potential is 
complex-valued. The agreement between the analytic and numerical results can be improved 
by using finer lattice spacings and larger number of points. 
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